
A dark morph of a male Fowler’s Toad (Anaxyrus fowled) in a cranberry (Vaccinium macrocarpon ) bog. The toad was vocalizing 
just before this picture was taken. The uninflated and darkly pigmented vocal sac is partially visible. Photo by Brad Timm. 
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Abstract .—We assessed the effects of landscape structure on Anaxyrus fowleri (Fowler’s Toad) site 
occupancy using 14 years of call survey data collected from 250 sites in Virginia and Maryland, and 
landscape variables derived from the National Wetlands Inventory, U.S. Census Bureau, National 
Land Cover Databases, and U.S. Department of Agriculture. We also conducted a time series 
analysis on A. fowleri occupancy rates using call survey data collected throughout Virginia and 
Maryland. We found A. fowleri site occupancy to be negatively affected by deciduous forest, hay 
crops, development and agricultural pesticides, and we identified a negative interannual trend in 
occupancy rates between 1999 and 2012. 
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Introduction 

Amphibian populations are declining globally, with an¬ 
thropogenic degradation of landscapes near wetlands 
having a major impact on many species of pond-breeding 
amphibians (Findlay and Houlahan 1997; Blaustein and 
Kiesecker 2002). Upland habitats surrounding wetlands 
are vital for successful dispersal, foraging, and non¬ 
breeding activities, making upland habitat quality critical 
to the life history of pond-breeding amphibians (Wind- 
miller 1996; Semlitsch 2000; Gibbons 2003; Bartlet et 
al. 2004). The negative effects of landscape degradation 
on many amphibian species, including hydroperiod al¬ 
teration, pollution of wetlands from roadway runoff and 
agricultural chemicals, and mechanical disturbance of 
foraging, retreat and burrowing sites (Luo et al. 1999; 
Turtle 2001; Gray et al. 2004b) are fairly well under¬ 
stood. These anthropogenic disturbances can ultimately 
impact mobility and survival of larval, juvenile, and adult 
amphibians, and can lead to population declines and ex¬ 
tirpations (Blaustein and Kiesecker 2002; Gray et al. 
2004a). Thus, understanding how natural and anthropo- 
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genic landscape-level processes effect amphibian popu¬ 
lations is critical to amphibian conservation. 

The Fowler’s Toad, Anaxyrus (Bufo) fowleri, is wide¬ 
ly but irregularly distributed throughout the eastern Unit¬ 
ed States, occurring from southern New England to the 
Florida Panhandle and as far west as Missouri, Arkansas, 
and Louisiana (Netting and Goin 1945; Green 1992; Kl- 
emens 1993; Conant and Collins 1998). Though typically 
associated with coastal dune systems and scrub-pine for¬ 
ests, A. fowleri also occurs in rocky and sparsely veg¬ 
etated areas in dry, sandy, deciduous woodlands, and in 
agricultural and developed areas (Schlaugh 1976, 1978; 
Klemens 1993; Zampella and Bunnell 2000; Rubbo and 
Kiesecker 2005; Gooch et al. 2006). Some biologists and 
naturalists have thus described A. fowleri as being toler¬ 
ant of urbanization, and scarification in agricultural ar¬ 
eas (Ferguson 1960; Martof et al. 1980; Klemens 1993; 
Rubbo and Kiesecker 2005; Gooch et al. 2006). How¬ 
ever, other studies suggest that A. fowleri are habitat spe¬ 
cialists sensitive to environmental perturbations (B reden 
1988; Green 2005; Tupper and Cook 2008). 
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In Canada, A. fowleri are federally protected (Oldham 
2003), but they are not considered a species of concern in 
the United States. Anaxyrus fowleri populations are be¬ 
lieved to be relatively stable and abundant in the eastern 
United States (Conant and Collins 1998). However, A. 
fowleri extirpations have been documented in the north¬ 
eastern and southeastern United States (Breden 1988; 
Klemens 1993; Mierzwa et al. 1998; Tupper and Cook 
2008; Walls et al., 2011; Milko 2012). These extirpations 
were largely attributed to anthropogenic disturbances, 
such as habitat degradation, pesticide application, road 
mortality (National Park Service, unpubl. data), hydro¬ 
period alteration, competition from invasive species, and 
probable increased predation pressures from urban toler¬ 
ant predators such as skunks {Mephitis mephitis ) and rac¬ 
coons (Procyon lotor) (Schaff and Garton 1970; Lazell 
1976; Groves 1980; Klemens 1993; Tupper and Cook 
2008; Milko 2012). 

In the southern mid-Atlantic region, A. fowleri occur 
throughout Virginia and Maryland, but are less common 
outside of the Coastal Plain (Mitchell and Reay 1999). 
Coastal regions are thought to contain more favorable 
upland habitats for this species (see Martof et al. 1980; 
Mitchell and Reay 1999; Cook et al. In prep), but much 
of the Coastal Plain in Virginia and Maryland is more 
densely populated and intensely developed than west¬ 
ern regions. For instance, the mid-Atlantic Coastal Plain 
has the highest population density and second-highest 
growth rate of all ecoregions in Virginia (VGDIF 2005). 
If A. fowleri are sensitive to landscape perturbations, hu¬ 
man population growth and development may lead to A. 
fowleri population declines in the southern mid-Atlantic. 

To the best of our knowledge, quantitative data de¬ 
scribing the effects of landscape-level variables on A. 
fowleri are non-existent for the mid-Atlantic and are lim¬ 
ited elsewhere (see Gooch et al. 2006; Tupper and Cook 
2008; Birx-Raybuck 2010; Eskew et al. 2011). Occupan¬ 
cy trend analyses of A. fowleri populations indicate that 


they are stable in most mid-Atlantic states (except Mary¬ 
land; see Weir et al. 2014), but these analyses are tempo¬ 
rally limited (Weir et al. 2009). Thus, critical thresholds in 
landscape-level variables essential to A. fowleri occupan¬ 
cy are unknown and it is unclear if southern mid-Atlantic 
populations are stable over the long term. Therefore our 
objectives were to identify and describe landscape-level 
variables that influence A. fowleri site occupancy and to 
complete a more comprehensive time-series analysis for 
this species in the southern mid-Atlantic. 

Materials and Methods 

Site selection and data collection 

We randomly selected 250 sites in Virginia and Mary¬ 
land for landscape-level analyses (Fig. 1). Selected sites 
were North American Amphibian Monitoring Program 
(NAAMP) calling anuran survey points (adjacent to 
wetlands) that were surveyed with anuran call counts be¬ 
tween 1999 and 2012 (Weir and Mossman 2005). Move¬ 
ment data for A. fowleri are limited, but available studies 
indicate that a 1 km buffer surrounding breeding wetlands 
is a biologically meaningful distance for analyzing the 
effects of landscape features on anuran (including Anaxy¬ 
rus spp.) occurrence (Clarke 1974; Miaud et al. 2000; 
Muths 2003; Bartlet et al. 2004; Smith and Green 2005; 
Forester et al. 2006). Therefore, our landscape variables 
were derived from 1 km buffers surrounding calling sur¬ 
vey points. Any calling survey points found to have over¬ 
lapping buffers were removed from analysis. Anuran call 
data (ranked ordinal values based on chorus intensity [0- 
3]) and sampling covariates (ambient temperatures, sky 
and wind conditions, and noise disturbance levels) were 
collected in accordance with NAAMP guidelines by 
trained NAAMP volunteers (Weir and Mossman, 2005). 



Fig. 1 . Map of calling anuran surveys conducted in Maryland and Virginia. Closed circles indicate sites occupied by A. fowleri and 
open circles indicate unoccupied sites. 
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Calculation of landscape variables 

We quantified landscape variables using data from four 
publicly available sources: (1) the National Wetlands 
Inventory (NWI) from the U.S. Fish and Wildlife Ser¬ 
vice; (2) 2012 TIGER/Line road data from the U.S. 
Census Bureau; (3) the National Land Cover Database 
(NLCD2006) from the Multi-Resolution Land Charac¬ 
teristics Consortium; and (4) the National Pesticide Use 
Database from the U.S. Department of Agriculture. Ini¬ 
tial data manipulation was done in QuantumGIS (QGIS; 
QGIS Development Team 2011). For NWI data, we ex¬ 
tracted distance from calling survey sites to nearest wet¬ 
land and determined the number and types of wetlands 
within a 1 km buffer zone of calling survey sites. Using 
TIGER/Line road files, we calculated road length and 
type within 1 km buffers. We prepared land cover data 
by clipping NLCD2006 data for each buffer into an indi¬ 
vidual raster file. These files were then imported into R 
(R Core Team 2013) and analyzed using the SDMTools 
package (VanDerWal 2013). Total pesticide application 
rate (kg/km 2 ) at each site was determined by calculating 
the sum of application levels within each buffer for all 
pesticides listed in the National Pesticide Use Database. 

Data Analyses 

We used the R package Unmarked (Fiske and Chandler 
2011) to identify landscape-level variables associated 
with A. fowled occupancy. Landscape level habitat data 
were only available for a single year, so we used MacK¬ 


enzie et al.’s (2002) occupancy model to account for 
imperfect detectability, particularly false-negative detec¬ 
tions, when evaluating habitat variables. False-positive 
detections can also result in high site occupancy biases, 
but false-positive detection rates vary between species, 
and previous tests of NAAMP volunteers resulted in no 
false-positive detections of A. fowled , even amongst 
inexperienced volunteers (Genet and Sargent 2003; Ro- 
yle and Link 2006; McClintock et al. 2010). Therefore, 
false-positive detections of A. fowled were unlikely to be 
at levels high enough to bias site occupancy and were not 
considered in the modeling process. 

We assessed models using a multimodel inference 
approach (see Burnham et al. 2011). We used 19 non- 
correlated site covariates (Table 1) considered to be bio¬ 
logically meaningful in anuran breeding site selection 
when creating a priori models (Cushman 2006). Julian 
date (date) and temperature (temp) were found to af¬ 
fect detection probability; therefore these two sampling 
covariates were used in all models. We ranked compet¬ 
ing models with Akaike Information Criterion (for data 
sets with high independent to dependent variable ratio 
[AICc]) by calculating differences between candidate 
models and the lowest AICc (A. AICc) model. We used 
Akaike weight (w.) for each model to guide selection. 

To determine change in occupancy between years, 
we fit a colonization-extinction model (MacKenzie et al. 
2003) using date and temp as covariates to account for 
differences in repeated sampling periods. We then used 
a smoothed trajectory to determine mean occupancy for 
each year (Weir et al. 2009). Serial autocorrelation in 


Table 1 . Landscape variables used in analyses of A. fowled calling anuran survey data for in Virginia and Maryland. 

Variable Description 


Crops 

Dec 

Dev 

Ever 

For 

Grass 

H 

Hay 

Mix 

Patch 

Pesticides 

Road 

Shrub 

Wavg 

Wdis 

Wet 

Wnear 

Wnum 

Wtype 


Proportion of area used for annual crops or perennial woody crops 
Proportion of forest with >75% canopy cover of deciduous trees 
Proportion of area that has been developed, including suburban and urban areas 
Proportion of forest with >75% or more evergreen trees 
Proportion of all forest habitats (Dec + Ever + Mix) 

Proportion of area with graminoids or herbaceous vegetation covering over 80% of land which might be grazed but not tilled 
Habitat Diversity, calculated as Shannon’s diversity index using habitat proportions 
Proportion of area planted with grass/legume mixtures used for grazing or hay crops 
Proportion of mixed forest with neither deciduous nor evergreen dominant 

Number of terrestrial habitat patches divided by total number of possible habitat patches (i.e., if each raster square represented a 
different type of habitat) 

Total kg/km 2 of agricultural pesticides applied within 1 km radius of buffer from CSS 
Total length of all roads in a 1 km radius buffer from CSS 

Proportion of area with canopy less than 5 m tall (e.g., shrubs and early successional forest) 

Average size of wetlands in a 1 km radius buffer from CSS ± in m 2 
Distance (m) of nearest wetland from Calling Survey Site (CSS) 

Proportion of total buffer area covered by wetlands 

Size (ha) of the wetland nearest to the CSS 

Number of wetlands in a 1 km radius buffer from CSS 

Number of different types of wetlands in a 1 km radius buffer from CSS 
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the residuals violated assumptions of a parametric lin¬ 
ear regression analysis; therefore we used a non-seasonal 
Autoregressive Integrated Moving Average (ARIMA) 
analysis to better understand changes in A. fowled oc¬ 
cupancy rates over time. We determined significance of 
parameters using a conditional least squares estimation. 
We assessed model fit with a Ljung-Box Q-test whereby 
a high P-value indicates that autocorrelation functions 
are not significantly different than white noise (Ljung 
and Box 1978). 

We used an empirical Bayesian approach to deter¬ 
mine conditional distribution of occurrence from the 
colonization-extinction model and then extrapolated best 
unbiased predictions of occupancy probability at each 
site. All sites with an occupancy probability > 0.75 were 
considered occupied. To determine distributions of spe¬ 
cific variables at occupied and unoccupied sites, we cre¬ 
ated site occupancy accumulation curves for each habitat 
variable. Habitat recommendations are based on maxi¬ 
mum values found at 90% of occupied sites. 

Trend analysis was completed in Minitab v.16 (Minit¬ 
ab Inc., Pennsylvania, USA) and all other analyses were 
completed in R v.3.0.2. Maps and figures were created 
using QGIS 2.0 and Excel 2010 (Microsoft Corp., Wash¬ 
ington, USA). 

Results 

Landscape analyses 

Two hundred fifty sites were sampled between 1999 and 
2012 throughout Maryland and Virginia, 108 (42.8%) of 
which had at least one detection of A. fowled. Approxi¬ 
mately eight-percent (300/3841) of sampling events re¬ 
sulted in detections. Not all sites were sampled every 
year (ranging from 1-14 years, x = 6.22; + 0.197), but 
most sites were surveyed > 10 times (Jr = 23.1; ± 1.09). 


Table 2. Top five models and full model from unmarked analy¬ 
sis of A. fowled occupancy in Virginia and Maryland. Mod¬ 
els are ranked from lowest to highest with AICc values. Julian 
date and temperature were used as sampling covariates in all 
models. The full model was constructed using the maximum 
number of site covariates which would create a model that con¬ 
verged: Wnum + Wavg + Road + Ever + Dec + Mix + Crops + 
Hay + Dev + Wet + Core + Patch + Pesticides. 


Model 

AICc 

A, AICc 

w i 

Dec + Hay + Dev + Pesticides 

1768.65 

0 

0.62 

Dec + Hay + Dev 

1769.55 

0.91 

0.38 

Dec + Hay + Pesticides 

1780.54 

11.90 

0.00 

Dec + Dev + Pesticides 

1791.16 

22.52 

0.00 

Pesticides 

1964.45 

195.80 

0.00 

Full model 

2958.68 

1190.03 

0.00 


Table 3. Transformed and untransformed beta coefficients from 
top occupancy models. 


Variable 

Untransformed 

Transformed 

Estimated SE 

(Intercept) 

3.887 

0.980 

0.017 

Dec 

-6.148 

0.002 

0.002 

Hay 

-5.373 

0.005 

0.005 

Dev 

-4.814 

0.008 

0.009 

Pesticides 

-0.005 

0.499 

0.013 


Bayesian analysis indicated that only two sites where no 
detections occurred had an occupancy probability > 0.75. 

We found strong support for a model indicating that 
deciduous forest, hay, development and pesticides nega¬ 
tively influenced A. fowled occupancy (AICc = 1768.65, 
w. = 0.62; Tables 2 and 3). The buffers ranged from 
0-97% deciduous forest (x = 0.32; + 0.015), 0-74% 
hay (x = 0.18; + 0.011), 0-83% development (x = 0.12; 
± 0.009), and 0-535 kg pesticides applied (x = 55.7; + 
5.231). Ninety-percent of occupied sites were covered by 
less than 25% development, 35% hay, and 50% decidu- 
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Fig. 2. Accumulation curve showing the maximum proportions 
of habitat variables found at each occupied site. Ninety-percent 
of occupied sites were covered by less than 0.25 development, 
0.35 hay, and 0.50 deciduous forest cover. 


Fig. 3. Accumulation curve showing the maximum agricultural 
pesticide application rates at occupied sites. The maximum 
amount of pesticides applied at 90% of occupied sites was 165 
kg/km 2 . 
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ous forest (Fig. 2). The maximum amount of pesticides 
applied at 90% of occupied sites was 165 kg/km 2 (Fig. 3). 

Time-series analyses 

Using ARIMA analysis, we found that a single autore¬ 
gressive term was contributing to interannual changes 
in A. fowleri occupancy rates (t = 4.32, P < 0.001). We 
confirmed that the model was valid, with uncorrelated 
residuals (Q = 13.2, df = 9, P = 0.152). Trend analysis in¬ 
dicated a downward trend in occupancy rates over time, 
with occupancy decreasing from 55.3% in 1999 to 29.5% 
in 2012 (Fig. 4). 


L 



£ 0.2 


U - 

IW 2000 2002 2004 2006 2flOK 2010 2012 

Year 


Discussion 
Forest cover 


Fig. 4. Mean site occupancy and confidence intervals of A. 
fowleri in Virginia and Maryland from 1999-2012. Time se¬ 
ries analysis indicates a 53% decrease in site occupancy, from 
55.3% in 1999 to 29.5% in 2012. 


Many contemporary landscape-level studies indicate 
a positive relationship between amphibian species and 
forested habitat within buffers of varying sizes around 
breeding ponds (see Cushman 2006). We were able to 
distinguish between forest types on a fairly large scale 
and identified a negative relationship between A. fowleri 
occupancy and deciduous forest. Although both species 
can be sympatric, A. fowleri is largely replaced by Amer¬ 
ican Toads (,Anaxyrus americanus ) in later successional 
forests that are dominated by moister, more nutrient rich 
soils, and hardwood trees (Wright and Wright 1967; 
Lazell 1976; Klemens 1993). Our results confirm long¬ 


standing observations made across A. fowlerf s range 
that suggest they are more common in early successional 
habitats that are either relatively open or dominated by 
mixed or coniferous forest (Hubbs 1918; Hoopes 1930; 
Netting and Goin 1945; Littleford 1946; Cory and Man- 
ion 1955; Wright and Wright 1967; Clarke 1974; Green 
1989; Lazell 1976; Klemens 1993; Zampella and Bun¬ 
nell 2000; Tupper and Cook 2008). 

Hay 

The proportion of area covered by grass/legume mixtures 
used for grazing or hay crops within the 1 km buffer was 


—■^ 



The Provincelands of Cape Cod National Seashore, Barnstable County, Massachusetts, USA. The reddish vegetation in the center 
of the photo is a cranberry ( Vaccinium macrocarpon ) bog, a wetland used for breeding by the Fowler’s toad. The surrounding land¬ 
scape is ideal for the Fowler’s toad and supports one of the largest populations of this species in the United States. The landscape 
contains a patchwork of sand, pitch pine ( Pinus rigida ), scrub oak ( Quercus ilicifolia ), and dune grass (Ammophila breviligulata). 
Photo by Rebecca Flaherty. 
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found to have a significantly negative impact on A. fowl¬ 
ed occurrence. Agricultural development can negatively 
affect anuran dispersal abilities, and soil compaction as¬ 
sociated with agricultural landscape alterations may pro¬ 
hibit anuran burrowing (Whalley et al. 1995; Jansen et 
al. 2001; Gray et al. 2004b). Wetlands within agricultural 
landscapes may be altered physically and biologically 
such that postmetamorphic anurans emerge smaller and 
presumably less fit (Beja and Alcazar 2003; Gray et al. 
2004a,b). ’ 

Development 

We found development to be another significant variable 
negatively affecting A. fowled occurrence. Development 
contributes to reduced genetic diversity in pond breeding 
amphibians, increased pollution of upland and wetland 
habitats, increased road mortality, and microclimate al¬ 
teration of remaining habitat patches (Soule 1987; Reh 
and Seitz 1989; Falirig 1995; deMayndier and Hunter 
1998; deMaynadier and Hunter 2000; Turtle 2001; Timm 
and McGarigal 2014). Various studies indicate that de¬ 
velopment and fragmentation is detrimental to amphib¬ 
ian persistence (see review in Cushman 2006), including 
A. fowled and congener A. amedcanus (Schlauch 1976, 
1978; Gibbs et al. 2005; Walls et al. 2011). However, 
studies conducted throughout the United States (e.g. 
New Jersey [Zampella and Bunnell 2000J, Pennsylvania 
[Rubbo and Kiesecker 20051, North Carolina (Gooch 
et al. 2006], and Louisiana [Milko 2012]) suggest that 
A. fowled are urban tolerant. While these studies have 
shown that A. fowled can occur in developed habitats, 
they do not include pre-development population sizes, 
are temporally and spatially limited, and are likely refer¬ 
ring to suburbanization rather than large-scale urbaniza¬ 
tion (see Schlauch 1978). 

Pesticides 

While the differences in AICc values between the top 
two models were small, four out of the top five mod¬ 
els included agricultural pesticide application levels as 
a negative covariate, indicating that pesticide exposure 
may play an important role in A. fowled site occupancy. 
Data suggest that chemical pesticides associated with ag¬ 
riculture have wide-ranging direct effects on amphibians, 
including endocrine disruption, immunosuppression, de¬ 
velopmental delays, and increased mortality (Mann et al. 
2009). Exposure to insecticides was found to be highly 
toxic to larval A. fowled in laboratory studies (see review 
in Green 2005) and had sub-lethal effects its congener, A. 
amedcanus, causing eye and limb deformities, increased 
time to metamorphosis, and reduced post-metamorphic 
body size (Harris et al. 2000; Boone and James 2003; 
Howe et al. 2004). Agricultural runoff containing pes¬ 
ticides may be contaminating certain wetlands in this 


study, thus potentially accounting for reduced A. fowled 
occupancy rates in agricultural landscapes. 

Trends 

By pooling data from Maryland and Virginia, we esti¬ 
mate that A. fowled occupancy has decreased by approx¬ 
imately 53% over the last 14 years. Wen- et al. (2009) 
found a significant, but negligible, negative occupancy 
trend for A. fowled in Delaware, and indicated unchang¬ 
ing occupancy rates in Virginia, Maryland, West Virgin¬ 
ia, and New Jersey between 2001 and 2007. However, a 
more recent study (conducted as the same time as ours, 
with a similar data set, see Weir et al. 2014) also indicat¬ 
ed A. fowled declines in Maryland. Differences in trend 
estimates between our study and Weir et al. (2009, 2014) 
may be due to differing sample sizes. Since we nearly 
doubled the scope of analysis of Weir et al. (2009) and 
have three more years than Weir et al. (2014), we be¬ 
lieve our results more accurately describe trends in A. 
fowled occupancy in Maryland and Virginia. Although 
a more comprehensive analysis is needed to identify the 
proximate causes of decline in A. fowled in Maryland 
and Virginia, we suspect that its declines are in part due 
to the recent loss of subclimax communities. Virginia 
has lost 51.6% of its softwood forest since 1940 (VDOF 
2014) and early successional habitats have been steadily 
juxtaposing to later successional seres throughout Mary¬ 
land: As of 2008 less than 10% of existing Maryland for¬ 
ests were occupied by early successional regimes (Lister 
2011 ). 

Anaxyrus fowled may be able to persist longer in mod¬ 
erately developed coastal environments (Schlauch 1978) 
than other pond-breeding amphibians due to their high 
fecundity rates, salt tolerance, desiccation resistance, 
and ability to breed in wetlands with varying hydrope¬ 
riod regimes (Wright and Wright 1967; Claussen 1974; 
Markow 1997; Tupper and Cook 2008; Birx-Raybuck 
2010; Eskew et al. 2012). The ability of A. fowled to oc¬ 
cupy these types of habitats is advantageous because they 
harbor fewer interspecific amphibian competitors (see 
Martof et al. 1980; Klemens 1993; Mitchell and Reay 
1999). If early successional habitats continue to become 
less widespread in the mid-Atlantic, coastal regions may 
become more important to the long-term persistence of 
A. fowled. However, much of the southern mid-Atlantic 
coastal plain is densely populated (VGDIF 2005) and in¬ 
tensely developed. Thus, successional changes occurring 
further inland coupled with increased urbanization of the 
southern mid-Atlantic coastal plain could potentially ex¬ 
acerbate declines. 

Conclusions 

Amphibian populations are more vulnerable to habitat 
loss and fragmentation when located on the margins of 
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their geographic range (Swihart et al. 2003). Our data in¬ 
dicate that even in the middle of their range, A. fowleri 
occupancy rates are declining. Landscapes most appro¬ 
priate for this species appear to contain only moderate 
amounts of deciduous forest (< 50%), few hay crops (< 
35%), relatively little development (< 25%), and low 
pesticide application rates. 
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